Подготовка к работе

Нам понадобятся следующие пакеты.

2. Подготовим данные

Проверим отсутсвующие значения и выбросим их при необходимости.

numberOfNA <- length(which(is.na(Boston)==T))
if(numberOfNA>0) {
  Boston <- Boston[complete.cases(Boston),]
}

Создадим данные для тестирования модели.

set.seed(123)
split <- sample.split(Boston,SplitRatio = 0.75) #assigns booleans to a new coloumn based on the split ratio
train <- subset(Boston,split==TRUE)
test <- subset(Boston,split==FALSE)

3. Разведочный анализ.

Посмотрим структуру датасета и базовую статистику.

str(Boston)
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
summary(Boston)
##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

Здесь мы можем видеть, что переменные «crim» и «black» принимают широкий диапазон значений.

Переменные «crim», «zn», «rm» и «black» имеют большую разницу между их медианной и средним значением, что указывает на множество выбросов в соответствующих переменных.

par(mfrow = c(1, 4))
boxplot(Boston$crim, main='crim',col='Sky Blue')
boxplot(Boston$zn, main='zn',col='Sky Blue')
boxplot(Boston$rm, main='rm',col='Sky Blue')
boxplot(Boston$black, main='black',col='Sky Blue')

Как и было предсказано, мы видим много выбросов.

Посчитаем корреляцию между признаками.

corr_matrix <- cor(Boston)
corrplot.mixed(corr_matrix)

Мы видим крайне сильную корреляцию - 0,9 - между индексом доступности радиальных магистралей (rad) и полной ставкой налога на имущество (tax). Также мы можем видеть сильную отрицательную корреляцию между расстоянием до центров занятости (dis) и долей некоммерческого бизнеса (indus), концентрацией оксидов азота (nox), долей старых домов (age); и между более низким статусом населения (lstat) и средней стоимостью домов, занимаемых владельцами (medv). И сильная положительная корреляция между indus и nox, indus и tax, nox и age, rm и medv.

Также мы видим мультиколлинеарность некоторых предикторов, что может стать проблемой при построении модели, поскольку даже небольшие изменения в любом из предикторов могут вызвать большие изменения в моделях.

4. Построение модели, прогнозирующей стоимость дома

Посмотрим распределение переменной medv.

qplot(x=medv,data=Boston,geom='histogram')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

На данной гистограмме распределения средней цены дома мы видим, что данные немного смещены влево. Однако, линейную регрессию нельзя строить на ненормально распределенных данных. Чтобы исправить распределние прологарифмируем параметр medv.

qplot(x=log(medv),data=Boston,geom='histogram')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

После логарифмирования у нас остался небольшой сдвиг вправо, но распределение выглядит более “нормальным”.

Задание 1. Построение полной модели.

Стандартизируем предикторы:

# сделаем наличие реки фактором 
Boston$chas <- as.factor(Boston$chas)
stand_Boston <- Boston %>% mutate_at(c('crim', 'zn', 'indus', 'nox', 'age', 'dis', 'tax', 'ptratio', 'black', 'lstat'), scale)
model1_stand = lm(log(medv)~., data= stand_Boston)
summary(model1_stand)
## 
## Call:
## lm(formula = log(medv) ~ ., data = stand_Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73361 -0.09747 -0.01657  0.09629  0.86435 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.320438   0.104407  22.225  < 2e-16 ***
## crim        -0.088351   0.011315  -7.808 3.52e-14 ***
## zn           0.027345   0.012815   2.134 0.033349 *  
## indus        0.016923   0.016886   1.002 0.316755    
## chas1        0.100888   0.034486   2.925 0.003598 ** 
## nox         -0.090199   0.017717  -5.091 5.07e-07 ***
## rm           0.090833   0.016728   5.430 8.87e-08 ***
## age          0.005928   0.014883   0.398 0.690567    
## dis         -0.103364   0.016811  -6.149 1.62e-09 ***
## rad          0.014267   0.002656   5.373 1.20e-07 ***
## tax         -0.105466   0.025368  -4.157 3.80e-05 ***
## ptratio     -0.082856   0.011337  -7.309 1.10e-12 ***
## black        0.037757   0.009815   3.847 0.000135 ***
## lstat       -0.207344   0.014496 -14.304  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1899 on 492 degrees of freedom
## Multiple R-squared:  0.7896, Adjusted R-squared:  0.7841 
## F-statistic: 142.1 on 13 and 492 DF,  p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(model1_stand)

Полная модель:

  • R-squared составляет 0.7841

  • F-statistic составляет 142.1

Полная модель объясняет 78% изменчивости. Протестируем, есть ли параметры, которые незначительно влияют на модель. Для этого воспользуемся критерием Акаике.

step = stepAIC(model1_stand,direction = 'both')
## Start:  AIC=-1667.19
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## - age      1    0.0057 17.755 -1669.0
## - indus    1    0.0362 17.786 -1668.2
## <none>                 17.749 -1667.2
## - zn       1    0.1643 17.914 -1664.5
## - chas     1    0.3088 18.058 -1660.5
## - black    1    0.5339 18.283 -1654.2
## - tax      1    0.6235 18.373 -1651.7
## - nox      1    0.9351 18.684 -1643.2
## - rad      1    1.0413 18.791 -1640.3
## - rm       1    1.0637 18.813 -1639.7
## - dis      1    1.3639 19.113 -1631.7
## - ptratio  1    1.9270 19.676 -1617.0
## - crim     1    2.1995 19.949 -1610.1
## - lstat    1    7.3809 25.130 -1493.2
## 
## Step:  AIC=-1669.03
## log(medv) ~ crim + zn + indus + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## - indus    1    0.0363 17.791 -1670.0
## <none>                 17.755 -1669.0
## + age      1    0.0057 17.749 -1667.2
## - zn       1    0.1593 17.914 -1666.5
## - chas     1    0.3138 18.069 -1662.2
## - black    1    0.5431 18.298 -1655.8
## - tax      1    0.6205 18.376 -1653.7
## - nox      1    0.9645 18.720 -1644.3
## - rad      1    1.0356 18.791 -1642.3
## - rm       1    1.1452 18.900 -1639.4
## - dis      1    1.5471 19.302 -1628.8
## - ptratio  1    1.9224 19.677 -1619.0
## - crim     1    2.1988 19.954 -1612.0
## - lstat    1    8.1949 25.950 -1479.0
## 
## Step:  AIC=-1670
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## <none>                 17.791 -1670.0
## + indus    1    0.0363 17.755 -1669.0
## + age      1    0.0058 17.786 -1668.2
## - zn       1    0.1451 17.936 -1667.9
## - chas     1    0.3399 18.131 -1662.4
## - black    1    0.5344 18.326 -1657.0
## - tax      1    0.6139 18.405 -1654.8
## - nox      1    0.9350 18.726 -1646.1
## - rad      1    1.0088 18.800 -1644.1
## - rm       1    1.1171 18.909 -1641.2
## - dis      1    1.7385 19.530 -1624.8
## - ptratio  1    1.8862 19.678 -1621.0
## - crim     1    2.2229 20.014 -1612.4
## - lstat    1    8.1604 25.952 -1481.0
step$anova
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat
## 
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
## 
##      Step Df    Deviance Resid. Df Resid. Dev       AIC
## 1                              492   17.74938 -1667.194
## 2   - age  1 0.005723781       493   17.75510 -1669.031
## 3 - indus  1 0.036264380       494   17.79137 -1669.999
step_model = lm(log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat,data=stand_Boston)

summary(step_model)
## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     rad + tax + ptratio + black + lstat, data = stand_Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73400 -0.09460 -0.01771  0.09782  0.86290 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.328996   0.101070  23.043  < 2e-16 ***
## crim        -0.088757   0.011298  -7.856 2.49e-14 ***
## zn           0.025361   0.012637   2.007 0.045308 *  
## chas1        0.105148   0.034228   3.072 0.002244 ** 
## nox         -0.083634   0.016414  -5.095 4.97e-07 ***
## rm           0.090673   0.016281   5.569 4.20e-08 ***
## dis         -0.108878   0.015671  -6.948 1.18e-11 ***
## rad          0.013446   0.002540   5.293 1.82e-07 ***
## tax         -0.094023   0.022774  -4.129 4.28e-05 ***
## ptratio     -0.081025   0.011196  -7.237 1.77e-12 ***
## black        0.037679   0.009781   3.852 0.000133 ***
## lstat       -0.204262   0.013570 -15.053  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1898 on 494 degrees of freedom
## Multiple R-squared:  0.7891, Adjusted R-squared:  0.7844 
## F-statistic: 168.1 on 11 and 494 DF,  p-value: < 2.2e-16

Новая модель:

  • R-squared составляет 0.7864

  • F-statistic составляет 104.3

С помощью функции stepAIC мы последовательно рассчитали AIC для некоторых моделей и удалили предикторы, которые незначительно влияют на модель, что позволило несколько увеличить R-squared.

Задание 2. Проведем диагностику полной модели.

  1. Проверим на мультиколлениарность, которую мы заметили ранее.
library(car)
vif(step_model)
##     crim       zn     chas      nox       rm      dis      rad      tax 
## 1.789704 2.239229 1.059819 3.778011 1.834806 3.443420 6.861126 7.272386 
##  ptratio    black    lstat 
## 1.757681 1.341559 2.581984

Для дальнейшей работы нужно удалить предикторы с VIF > 5 (rad, tax).

  1. Проверим остатки:
plot(step_model$residuals)
abline(h=0,col='red')

Мы видим что остатки распределены ненормально.

  1. Найдем выбросы:
outlierTest(step_model)
##      rstudent unadjusted p-value Bonferroni p
## 413  4.775740         2.3647e-06    0.0011965
## 372  4.340845         1.7232e-05    0.0087195
## 373  4.007183         7.0941e-05    0.0358960
## 402 -3.952037         8.8814e-05    0.0449400

Мы видим 4 выброса, влияющие на нашу модель.

Наибольший вклад вносит переменная lstat - -0.207344.

Задание 3. Постройте график предсказаний стоимости от переменной, которая обладает

наибольшим по модулю коэффициентом.

mod_hi <- lm(medv ~ as.numeric(lstat), data = stand_Boston)

test$predicted.medv <- predict(mod_hi,test)
pl1 <-test %>% 
  ggplot(aes(medv,predicted.medv)) +
  geom_point(alpha=0.5) + 
  stat_smooth(aes(colour='black')) +
  xlab('Значения medv') +
  ylab('Предсказанные значения medv') +
  theme_bw()

ggplotly(pl1)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Дополнительное задание.

Продолжим преобразовывать модель на НЕстандартизованных данных.

model1_stand = lm(log(medv)~., data= stand_Boston)
summary(model1_stand)
## 
## Call:
## lm(formula = log(medv) ~ ., data = stand_Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73361 -0.09747 -0.01657  0.09629  0.86435 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.320438   0.104407  22.225  < 2e-16 ***
## crim        -0.088351   0.011315  -7.808 3.52e-14 ***
## zn           0.027345   0.012815   2.134 0.033349 *  
## indus        0.016923   0.016886   1.002 0.316755    
## chas1        0.100888   0.034486   2.925 0.003598 ** 
## nox         -0.090199   0.017717  -5.091 5.07e-07 ***
## rm           0.090833   0.016728   5.430 8.87e-08 ***
## age          0.005928   0.014883   0.398 0.690567    
## dis         -0.103364   0.016811  -6.149 1.62e-09 ***
## rad          0.014267   0.002656   5.373 1.20e-07 ***
## tax         -0.105466   0.025368  -4.157 3.80e-05 ***
## ptratio     -0.082856   0.011337  -7.309 1.10e-12 ***
## black        0.037757   0.009815   3.847 0.000135 ***
## lstat       -0.207344   0.014496 -14.304  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1899 on 492 degrees of freedom
## Multiple R-squared:  0.7896, Adjusted R-squared:  0.7841 
## F-statistic: 142.1 on 13 and 492 DF,  p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(model1_stand)

Новая модель:

  • R-squared составляет 0.7841

  • F-statistic составляет 142.1

step = stepAIC(model1_stand,direction = 'both')
## Start:  AIC=-1667.19
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## - age      1    0.0057 17.755 -1669.0
## - indus    1    0.0362 17.786 -1668.2
## <none>                 17.749 -1667.2
## - zn       1    0.1643 17.914 -1664.5
## - chas     1    0.3088 18.058 -1660.5
## - black    1    0.5339 18.283 -1654.2
## - tax      1    0.6235 18.373 -1651.7
## - nox      1    0.9351 18.684 -1643.2
## - rad      1    1.0413 18.791 -1640.3
## - rm       1    1.0637 18.813 -1639.7
## - dis      1    1.3639 19.113 -1631.7
## - ptratio  1    1.9270 19.676 -1617.0
## - crim     1    2.1995 19.949 -1610.1
## - lstat    1    7.3809 25.130 -1493.2
## 
## Step:  AIC=-1669.03
## log(medv) ~ crim + zn + indus + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## - indus    1    0.0363 17.791 -1670.0
## <none>                 17.755 -1669.0
## + age      1    0.0057 17.749 -1667.2
## - zn       1    0.1593 17.914 -1666.5
## - chas     1    0.3138 18.069 -1662.2
## - black    1    0.5431 18.298 -1655.8
## - tax      1    0.6205 18.376 -1653.7
## - nox      1    0.9645 18.720 -1644.3
## - rad      1    1.0356 18.791 -1642.3
## - rm       1    1.1452 18.900 -1639.4
## - dis      1    1.5471 19.302 -1628.8
## - ptratio  1    1.9224 19.677 -1619.0
## - crim     1    2.1988 19.954 -1612.0
## - lstat    1    8.1949 25.950 -1479.0
## 
## Step:  AIC=-1670
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## <none>                 17.791 -1670.0
## + indus    1    0.0363 17.755 -1669.0
## + age      1    0.0058 17.786 -1668.2
## - zn       1    0.1451 17.936 -1667.9
## - chas     1    0.3399 18.131 -1662.4
## - black    1    0.5344 18.326 -1657.0
## - tax      1    0.6139 18.405 -1654.8
## - nox      1    0.9350 18.726 -1646.1
## - rad      1    1.0088 18.800 -1644.1
## - rm       1    1.1171 18.909 -1641.2
## - dis      1    1.7385 19.530 -1624.8
## - ptratio  1    1.8862 19.678 -1621.0
## - crim     1    2.2229 20.014 -1612.4
## - lstat    1    8.1604 25.952 -1481.0
step$anova
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat
## 
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
## 
##      Step Df    Deviance Resid. Df Resid. Dev       AIC
## 1                              492   17.74938 -1667.194
## 2   - age  1 0.005723781       493   17.75510 -1669.031
## 3 - indus  1 0.036264380       494   17.79137 -1669.999
step_model = lm(log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat,data=stand_Boston)

summary(step_model)
## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     rad + tax + ptratio + black + lstat, data = stand_Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73400 -0.09460 -0.01771  0.09782  0.86290 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.328996   0.101070  23.043  < 2e-16 ***
## crim        -0.088757   0.011298  -7.856 2.49e-14 ***
## zn           0.025361   0.012637   2.007 0.045308 *  
## chas1        0.105148   0.034228   3.072 0.002244 ** 
## nox         -0.083634   0.016414  -5.095 4.97e-07 ***
## rm           0.090673   0.016281   5.569 4.20e-08 ***
## dis         -0.108878   0.015671  -6.948 1.18e-11 ***
## rad          0.013446   0.002540   5.293 1.82e-07 ***
## tax         -0.094023   0.022774  -4.129 4.28e-05 ***
## ptratio     -0.081025   0.011196  -7.237 1.77e-12 ***
## black        0.037679   0.009781   3.852 0.000133 ***
## lstat       -0.204262   0.013570 -15.053  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1898 on 494 degrees of freedom
## Multiple R-squared:  0.7891, Adjusted R-squared:  0.7844 
## F-statistic: 168.1 on 11 and 494 DF,  p-value: < 2.2e-16

Новая модель:

  • R-squared составляет 0.7844

  • F-statistic составляет 168.1

library(car)
vif(step_model)
##     crim       zn     chas      nox       rm      dis      rad      tax 
## 1.789704 2.239229 1.059819 3.778011 1.834806 3.443420 6.861126 7.272386 
##  ptratio    black    lstat 
## 1.757681 1.341559 2.581984
outlierTest(step_model)
##      rstudent unadjusted p-value Bonferroni p
## 413  4.775740         2.3647e-06    0.0011965
## 372  4.340845         1.7232e-05    0.0087195
## 373  4.007183         7.0941e-05    0.0358960
## 402 -3.952037         8.8814e-05    0.0449400

У нас есть 4 выброса, которые мы удалим.

Boston_1 = Boston[-c(413,372, 373, 402),,drop=T]# drop the points
row.names(Boston_1)=1:nrow(Boston_1)

model2 = lm(log(medv)~crim+zn+chas+nox+dis+ptratio+black+lstat,data=Boston)
summary(model2)
## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + dis + ptratio + 
##     black + lstat, data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70805 -0.12147 -0.01373  0.10405  0.82894 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7010244  0.1421502  33.071  < 2e-16 ***
## crim        -0.0080105  0.0012719  -6.298 6.65e-10 ***
## zn           0.0013826  0.0005565   2.484 0.013311 *  
## chas1        0.1302984  0.0363485   3.585 0.000371 ***
## nox         -0.7322562  0.1352610  -5.414 9.63e-08 ***
## dis         -0.0578514  0.0077684  -7.447 4.25e-13 ***
## ptratio     -0.0376285  0.0048687  -7.729 6.06e-14 ***
## black        0.0002958  0.0001120   2.642 0.008505 ** 
## lstat       -0.0353753  0.0017219 -20.545  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2025 on 497 degrees of freedom
## Multiple R-squared:  0.7584, Adjusted R-squared:  0.7545 
## F-statistic:   195 on 8 and 497 DF,  p-value: < 2.2e-16
plot(model2)

Новая модель:

  • R-squared составляет 0.7545
  • F-statistic составляет 195

Не смотря на то, что мы устранили мульитколлинеанрость и незначимые предикторы, предсказания модели стали немного хуже (75%).

vif(model2)
##     crim       zn     chas      nox      dis  ptratio    black    lstat 
## 1.473811 2.074443 1.049493 3.024847 3.294717 1.367971 1.286382 1.861558

Проведем диагностику модели

# Assumption - mean of the residuals is = 0 
mean(model2$residuals)
## [1] 4.05783e-18
plot(model2$residuals)
abline(h=0,col='red')

Остатки распределены более меннее нормально.

Найдем выбросы.

outlierTest(model2)
##     rstudent unadjusted p-value Bonferroni p
## 413 4.270130         2.3419e-05     0.011850
## 372 4.052312         5.8859e-05     0.029782

У нас есть два выброса, которые можно удалить.

Boston_1 = Boston[-c(413,372),,drop=T]# drop the points
row.names(Boston_1)=1:nrow(Boston_1)

model3 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
              ptratio + black + lstat,data=Boston_1)
par(mfrow=c(2,2))
plot(model3)

summary(model3)
## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     ptratio + black + lstat, data = Boston_1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70878 -0.09807 -0.00993  0.08974  0.77161 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.6741127  0.1896496  19.373  < 2e-16 ***
## crim        -0.0084826  0.0011748  -7.220 1.97e-12 ***
## zn           0.0007490  0.0005193   1.442 0.149874    
## chas1        0.1199066  0.0335829   3.570 0.000391 ***
## nox         -0.6147940  0.1258689  -4.884 1.40e-06 ***
## rm           0.1123588  0.0157923   7.115 3.96e-12 ***
## dis         -0.0447458  0.0073176  -6.115 1.97e-09 ***
## ptratio     -0.0333834  0.0045429  -7.348 8.36e-13 ***
## black        0.0004409  0.0001050   4.197 3.20e-05 ***
## lstat       -0.0289046  0.0018797 -15.377  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1869 on 494 degrees of freedom
## Multiple R-squared:  0.7937, Adjusted R-squared:  0.7899 
## F-statistic: 211.1 on 9 and 494 DF,  p-value: < 2.2e-16

Новая модель:

  • R-squared составляет 0.7899

  • F-statistic составляет 211.1

После удаления выбросов, мы видим, что значение R-squared увеличилось. Однако один из параметров перестал быть значимым.

model3_step=stepAIC(model3,Boston_1,direction = 'both')
## Start:  AIC=-1680.97
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black + 
##     lstat
## 
##           Df Sum of Sq    RSS     AIC
## <none>                 17.247 -1681.0
## - zn       1    0.0726 17.320 -1680.8
## - chas     1    0.4451 17.692 -1670.1
## - black    1    0.6151 17.862 -1665.3
## - nox      1    0.8329 18.080 -1659.2
## - dis      1    1.3054 18.552 -1646.2
## - rm       1    1.7673 19.014 -1633.8
## - crim     1    1.8202 19.067 -1632.4
## - ptratio  1    1.8853 19.132 -1630.7
## - lstat    1    8.2551 25.502 -1485.8
model3_step$anov
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black + 
##     lstat
## 
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black + 
##     lstat
## 
## 
##   Step Df Deviance Resid. Df Resid. Dev       AIC
## 1                        494   17.24699 -1680.969
summary(model3_step)
## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     ptratio + black + lstat, data = Boston_1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70878 -0.09807 -0.00993  0.08974  0.77161 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.6741127  0.1896496  19.373  < 2e-16 ***
## crim        -0.0084826  0.0011748  -7.220 1.97e-12 ***
## zn           0.0007490  0.0005193   1.442 0.149874    
## chas1        0.1199066  0.0335829   3.570 0.000391 ***
## nox         -0.6147940  0.1258689  -4.884 1.40e-06 ***
## rm           0.1123588  0.0157923   7.115 3.96e-12 ***
## dis         -0.0447458  0.0073176  -6.115 1.97e-09 ***
## ptratio     -0.0333834  0.0045429  -7.348 8.36e-13 ***
## black        0.0004409  0.0001050   4.197 3.20e-05 ***
## lstat       -0.0289046  0.0018797 -15.377  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1869 on 494 degrees of freedom
## Multiple R-squared:  0.7937, Adjusted R-squared:  0.7899 
## F-statistic: 211.1 on 9 and 494 DF,  p-value: < 2.2e-16

Новая модель:

  • R-squared составляет 0.7899

  • F-statistic составляет 211.1

par(mfrow=c(2,2))
plot(model3_step)

Мы видим, что нет разницы между model3 и model3_step, а также Р-квадрат имеет одинаковые значения. Теперь можем проанализировать выбросы модели.

outlierTest(model3)
##     rstudent unadjusted p-value Bonferroni p
## 369 4.303254         2.0307e-05     0.010235
## 372 4.266323         2.3833e-05     0.012012

Мы видим 2 выбросы, которые мы можем выбросить и протестировать модель еще раз.

Boston_2=Boston_1[-c(371,400,373),,drop=T]
row.names(Boston_2)=1:nrow(Boston_2)
model4 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
                 tax + ptratio + black + lstat,data=Boston_2)
summary(model4)
## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     tax + ptratio + black + lstat, data = Boston_2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70835 -0.10150 -0.00960  0.08772  0.79154 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.624e+00  1.899e-01  19.086  < 2e-16 ***
## crim        -7.880e-03  1.217e-03  -6.473 2.34e-10 ***
## zn           7.611e-04  5.218e-04   1.459  0.14530    
## chas1        1.062e-01  3.335e-02   3.185  0.00154 ** 
## nox         -5.614e-01  1.376e-01  -4.079 5.27e-05 ***
## rm           1.154e-01  1.550e-02   7.445 4.38e-13 ***
## dis         -4.340e-02  7.194e-03  -6.033 3.18e-09 ***
## tax         -6.026e-05  8.344e-05  -0.722  0.47055    
## ptratio     -3.255e-02  4.937e-03  -6.594 1.11e-10 ***
## black        4.396e-04  1.048e-04   4.195 3.23e-05 ***
## lstat       -2.858e-02  1.870e-03 -15.288  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1831 on 490 degrees of freedom
## Multiple R-squared:  0.7971, Adjusted R-squared:  0.7929 
## F-statistic: 192.4 on 10 and 490 DF,  p-value: < 2.2e-16

Новая модель:

  • R-squared составляет 0.7929

  • F-statistic составляет 192.4

par(mfrow=c(2,2))
plot(model4)

Мы видим, что мы снова увеличили мощность нашей модели, но потеряли значение по другой переменной. Мы можем запустить stepAIC и посмотреть, сможем ли мы его удалить.

outlierTest(model4)
##      rstudent unadjusted p-value Bonferroni p
## 369  4.527033         7.5217e-06    0.0037683
## 371  4.517479         7.8548e-06    0.0039352
## 372  3.982469         7.8571e-05    0.0393640
## 398 -3.953793         8.8281e-05    0.0442290

Мы видим 4 выбросы, которые мы можем выбросить и протестировать модель еще раз.

Boston_3 =Boston_2[-c(369,371,372,398),,drop=T]
row.names(Boston_3)=1:nrow(Boston_3)
Boston_3$chas <- as.factor(Boston_3$chas)
model5 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
                 tax + ptratio + black + lstat,data=Boston_3)
summary(model5)
## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     tax + ptratio + black + lstat, data = Boston_3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64945 -0.09613 -0.01079  0.08989  0.69041 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.421e+00  1.781e-01  19.207  < 2e-16 ***
## crim        -7.826e-03  1.132e-03  -6.915 1.48e-11 ***
## zn           5.723e-04  4.854e-04   1.179 0.239025    
## chas1        7.966e-02  3.144e-02   2.533 0.011610 *  
## nox         -5.005e-01  1.282e-01  -3.905 0.000108 ***
## rm           1.385e-01  1.468e-02   9.436  < 2e-16 ***
## dis         -3.810e-02  6.719e-03  -5.670 2.45e-08 ***
## tax         -1.143e-04  7.801e-05  -1.465 0.143473    
## ptratio     -3.225e-02  4.591e-03  -7.023 7.35e-12 ***
## black        4.324e-04  9.785e-05   4.418 1.23e-05 ***
## lstat       -2.677e-02  1.789e-03 -14.961  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1701 on 486 degrees of freedom
## Multiple R-squared:   0.82,  Adjusted R-squared:  0.8163 
## F-statistic: 221.4 on 10 and 486 DF,  p-value: < 2.2e-16

Новая модель:

  • R-squared составляет 0.8163

  • F-statistic составляет 221.4

par(mfrow=c(2,2))
plot(model5)

Мы улучшили еще немного нашу модель и теперь она объясняет 82% изменчивости. Нет предела совершенству, меня устраивает модель, объясняющая 82% изменчивости. Наибольший вклад вносят параметры: количество комнат, количество оксида азота, отношение учеников к учителям. Так что я бы посоветовала клиенту, строить большие дома рядом со школой и с парком неподалеку.

Финальный анализ модели.

residuals <- data.frame('Residuals' = model5$residuals)
res_hist <- ggplot(residuals, aes(x=Residuals)) + geom_histogram(color='black', fill='skyblue') + ggtitle('Histogram of Residuals')
res_hist
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Остатки распределены нормально.

plot(model5, col='Sky Blue')

# пересчитаем тестовый датасет
set.seed(123)
split <- sample.split(Boston,SplitRatio = 0.75) #assigns booleans to a new coloumn based on the split ratio
train <- subset(Boston,split==TRUE)
test <- subset(Boston,split==FALSE)

test$predicted.medv <- predict(model5,test)
pl1 <-test %>% 
  ggplot(aes(medv,predicted.medv)) +
  geom_point(alpha=0.5) + 
  stat_smooth(aes(colour='black')) +
  xlab('Значения medv') +
  ylab('Предсказанные значения medv') +
  theme_bw()

ggplotly(pl1)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'